---
title: "Simulation Code (for Final Submission)"
date: "`January 2024"
output: pdf_document
---

```{r, echo=FALSE}
library(ggplot2)
library(ggthemes)
library(gridExtra)
library(reshape2)
library(patchwork)
```


```{r,echo=FALSE}
K_original=rnorm(1000000)
S_original=rnorm(1000000)
alpha_k=0.6
alpha_s=0.5
step=0.05

K_optimal=K_original[K_original>= log(alpha_k/alpha_s)+ S_original]
S_optimal=S_original[K_original<  log(alpha_k+alpha_s)+ S_original]
K_quantile=quantile(K_original,probs=seq(0.05,0.95,step))
S_quantile=quantile(S_original,probs=seq(0.05,0.95,step))
K_optimal_quantile=quantile(K_optimal,probs=seq(0.05,0.95,step))
S_optimal_quantile=quantile(S_optimal,probs=seq(0.05,0.95,step))
quantile_data=data.frame(
  Rank=seq(0.05,0.95,step),
  Manager=K_optimal_quantile,
  Worker=S_optimal_quantile,
  Uncensored=K_quantile
)
SkillPlot=ggplot(melt(quantile_data,id="Rank"),aes(x=Rank,y=value,shape=variable))+
  geom_line()+geom_point()+ggtitle("Log Occupational Skills by Rank")+xlab("Rank")+ylab("Occupational Log Skill")+xlim(c(0,1))+theme_light()+theme(legend.position = "none")
```



```{r,echo=FALSE}
mylength=length(K_optimal_quantile)+1

pi=rep(NA,mylength)
w=rep(NA,mylength)

pi[1]=1/2
w[1]=1/2
for (K in 2:mylength){
  pi[K]=pi[K-1]+alpha_k* exp(alpha_k*K_optimal_quantile[K-1]+alpha_s*S_optimal_quantile[K-1])
}
for (S in 2:mylength){
  w[S]=w[S-1]+alpha_s*exp(alpha_k*K_optimal_quantile[S-1]+alpha_s*S_optimal_quantile[S-1])
}
log_pi=log(pi)
log_w=log(w)


Earnings_data=data.frame(
  Rank=seq(0,0.95,step),
  Manager=log_pi,
  Worker=log_w
)

EarningsPlot=ggplot(melt(Earnings_data,id="Rank"),aes(x=Rank,y=value,shape=variable))+
  geom_line()+geom_point()+ggtitle("Log Earnings by Rank")+xlab("Rank")+ylab("Log Earnings")+xlim(c(0,1))+theme_light()+theme(legend.position = "none")
```

```{r fig.height=6, fig.width=4}

SkillPlot_legend=ggplot(melt(quantile_data,id="Rank"),aes(x=Rank,y=value,shape=variable))+
  geom_line()+geom_point()+ggtitle("Skills by Rank")+xlab("Rank")+ylab("Occupational Log Skill")+xlim(c(0,1))+theme_light()+theme(legend.position = "bottom")+guides(shape = guide_legend(title = ""))

# Create user-defined function, which extracts legends from ggplots
extract_legend <- function(my_ggp) {
  step1 <- ggplot_gtable(ggplot_build(my_ggp))
  step2 <- which(sapply(step1$grobs, function(x) x$name) == "guide-box")
  step3 <- step1$grobs[[step2]]
  return(step3)
}

shared_legend <- extract_legend(SkillPlot_legend)

# Draw plots with shared legend
grid.arrange(arrangeGrob(SkillPlot, EarningsPlot, nrow = 2),
             shared_legend, nrow = 2, heights = c(10, 1))

```


